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Abstract. We give the equations governing the shear evolution in Bianchi 
spacetimes for general /(ij)-theories of gravity. We consider the case of i?"- 
gravity and perform a detailed analysis of the dynamics in Bianchi I cosmologies 
which exhibit local rotational symmetry. We find exact solutions and study their 
behaviour and stability in terms of the values of the parameter n. In particular, 
we found a set of cosmic histories in which the universe is initially isotropic, then 
develops shear anisotropics which approaches a constant value. 
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1. Introduction 

In the last few years there has been renewed interest in theories of gravity where 
the gravitational Lagrangian is a non-linear function of the scalar curvature. These 
f {R)-theories of gravity can take on a number of forms, the majority of the functions 
considered being of the type R + e_R™. Theories with rn = —1 have been proposed 
as possible alternatives to sources of dark energy to explain the observed cosmic 
acceleration [1,2]. Solar system experiments do however constrain these type of 
theories for any corrections higher than (quadratic gravity) [3]. In these theories 
corrections to the characteristic length scale of General Relativity (GR) are introduced 
through the addition of a new length scale which is determined by the constant e. 

There are however forms of f{R) which do not alter the characteristic length 
scale, for example i?", in which GR is recovered when n — 1. These i?"-gravity 
theories have many attractive features, such as simple exact solutions which allows 
for comparison with observations [4,5]. There are however some caveats, in particular 
the stability and global behaviour of the underlying cosmological model is not well 
understood. The dynamical systems approach [6] can address some of these problems, 
since it provides one with exact solutions through the determination of fixed points 
and a (qualitative) description of the global dynamics of the system. Carloni et al [7] 
have recently used this method to study the dynamics of _R"-theories in Friedmann- 
Lemaitre-Robertson- Walker (FLRW) universes. Clifton and Barrow [8] used the 
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dynamical systems approach to determine the extent to which exact solutions can 
be considered as attractors of spatially flat universes at late times. They compared 
the predictions of these results with a range of observations and were able to show 
that the parameter n in FLRW may only deviate from GR by a very small amount 
(n - 1 - 10-19). 

The dynamics of anisotropic models with /(i?)-gravity have not been studied 
as intensively as their FLRW counterparts and it is therefore not known how the 
behaviour of the shear is modified in these theories of gravity. Bianchi spacetimes 
with isotropic 3-surfaces have been investigated for the quadratic theory [9] and it 
was found that in Bianchi I cosmologies the universe isotropises slower than in the 
Einstein case. The equations governing the evolution of shear in Bianchi spacetimes 
for general /(i?)-theories can be found from the trace- free Gauss-Codazzi equations 
(see for example [10,11]). However, it is not easy to solve these equations since the 
shear depends non-linearly on the Ricci scalar. Consequently the dynamical systems 
approach provides us with the best means of understanding the dynamics of these 
models. 

In GR the vacuum Kasner solutions [12] and their fluid flUed counterparts, the 
Type I Bianchi models, proved useful as a starting point for the investigation of the 
structure of anisotropic models. Barrow and Clifton [13, 14] have recently shown that 
it is also possible to find solutions of the Kasner type for i?"-gravity models. In [13] 
they showed that exact Kasner-like solutions do exist in the range of parameter n 
for 1/2 < n < 5/4 but with different Kasner- index relations to the ones in GR. In 
this paper we extend the dynamical systems analysis of i?"-gravity [7] to Bianchi 
I cosmological models that exhibit local rotational symmetry (LRS) [15-17]. LRS 
spacetimes geometries are subgroups within anisotropic spacetimes in which isotropics 
can occur around a point within the spacetime in 1- or 3-dimensions. Thus there 
exists a unique preferred spatial direction at each point which constitutes a local axis 
of symmetry. All observations are identical under rotation about the axis and are the 
same in all spatial directions perpendicular to that direction [15, 16]. 

Our analysis of the vacuum case revealed some interesting features regarding 
the evolution of Bianchi I cosmologies in HOTG. The phase space contains one 
isotropic fixed point and a line of fixed points with shear. The isotropic fixed point 
is an attractor (stable node) for values of the parameter n in the ranges n < 1/2, 
1/2 < 71 < 1 and n > 5/4. In the range 1 < n < 5/4 this point is a repeller (unstable 
node) and therefore may be seen as a past attractor. The existence of an isotropic 
past attractor implies, as in the case of braneworld models [18-21], that we do not 
require special initial conditions for inflation to start since the cosmological singularity 
is FLRW. Another feature of these models is that the shear evolution is independent 
of the value of n and is the same for all fixed points on the line. 

The outline of this paper is as follows: In section 2 we give the basic equations 
for the kinematical and dynamical variables and we also state the field equations 
for general /(-R)-theories of gravity. In section 3 we convert these equations into an 
autonomous set of equations for the case of i?". We then proceed in section 4 and 5 
to analyse this system for models with vacuum and matter respectively. 

The following conventions will be used in this paper: the metric signature is 
( — h +-I-); Latin indices run from to 3; V represents the usual covariant derivatives 
which may be split (1-1- 3-covariantly) with the spatial covariant derivative being 
denoted by D and the time derivative by a dot; units are used in which c ~ SirG = 1. 
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2. Kinematics and dynamics of /(i?)-gravity 

2.1. Kinematical and dynamical quantities 

We first state the relevant equations governing relativistic fluid dynamics (see for 
e.g. [10, 11]). For any given fluid 4- velocity vector field the projection tensor 
hab = dab + MaMfe projects iuto the instantaneous rest-space of a comoving observer. 
The first covariant derivative can be decomposed as 

VaUb = - Ua Ufa + DaUb ^ - Ua Ub + ^ Q hab + CFab + ^ab , (1) 

where aab is the symmetric shear tensor [aab = '^(ab)^ <^abU^ ~ 0, a'^a = 0), uoab 
is the vorticity tensor {iOab = W[ah] i '^ab = 0) and iia is the acceleration vector 
{ua = u^VbUa). S is the volume expansion (0 = V au"") which defines a length scale 
a along the fiow lines via the standard relation Q = —. 

The matter energy-momentum tensor can be decomposed relative to w° in 
the form 

T^l = ^UaUb + qaUb+ Ua Qb + P hab + T^ab , (2) 

where is the relativistic energy density, p the isotropic pressure, q"" the energy 
flux (ciaU°' = 0) and i^ab the trace-free anisotropic pressure (tt^o = 0, i^ab = ttj^j), 
T^abU^ = 0), all relative to 

The conservation equations '^^T^b ~ ^7 can be split with respect to m° and hab- 

A + Daq" + e{fl +p)+ 2{Uaqn + KbTT^a) = 0, (3) 
qia) + Ijap ^ j^^^ab + 4 q ^ ^a^ ^ + p) il" - Ub Tt"'' + C^^^fc = 0. (4) 

2.2. f{R) -gravity Field Equations 

In the case where the gravitational Lagrangian is a non-linear function of the scalar 
curvature f{R), the action reads 

A = J dx^^f{R)+ J CAidx^, (5) 

where Cm is the Lagrangian of the matter fields. The fourth order field equations can 
be obtained by varying 

T„Y = f'Rab - \fgab + {f'gab - glg'b) , (6) 

where primes denote derivatives with respect to R and Sab = VaVf,/'(i?). Sab can be 
decomposed as 

Sab = f'DaDbR + f'DaR DbR - f"'R (DbR Ua + DaR Ub) + fie UaUb 

-f" [OaR Ub + U^VciDbR) Ma] + /" [R UaUb ~ R{DaUb - UaUb)\ , (7) 

and so the D'Lambertian can be given by 

S = □/'(i?) = 5''^VaVfc/'(i?) 
= fD'-D^R + /'"D^R DcR + fu'D^R - /" {R + RQ)- f"'R^ . (8) 
The field equation (jHJ can be rewritten in the standard form 

Rab - \gabR = T'^b''^ , (9) 
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(when f'{R) 7^ 0) where the effective stress energy momentum tensor T^^^ is given 

by 

= [T^l + hab if f'R) + Sea {9^9% - 5'=^^)] ■ (10) 

As pointed out in [7] , no matter how compUcated the effective stress energy momentum 
tensor for the HOTG system is, it is always divergence free if V''Tj^j^ = 0. The total 
conservation equations therefore have the same form as those for standard matter and 
can thus be represented by © and Q). 

The higher order field equations may then be split (see [9, 10, 22]) to give the 
following contributions: 

= /'"M3p-/i + 2/-35], (11) 
Rabu'^u' = f [m -y + h'^'Sab] , (12) 

Rabu'^hl = /' -1 [-q, + Sabu'^h':] , (13) 

Rabh^X = /' Kd -{p+y + S) + Sabh'^M ■ (14) 

The propagation and constraint equations for /(i?)-theories are given by Ripple 
et al. [22]. In what follows we will use the Raychaudhuri equation 

e - Dau"" + 5 - (w,u«) + 2a^ -2^"" + f'-^[ii-\f + h''''Sab] = 0, (15) 

and the trace free Gauss-Codazzi equation, which holds for an irrotational matter fluid 
flow {uJab = 0): 

^Rab - ^CR)hab = - ^{ab) - 'd O^ab + D i^^Ub) + ^i^a^b) + f '^T^ab 

+ /' -1 [hlht - Ihabh'^"] 5ed, (16) 

where the 3-Ricci scalar is given by 

^ - + /' -1 + 3p + / - ^ + 2h''-'S,d\ ■ (17) 

3. Shear dynamics in Bianchi I cosmologies 

3.1. f{R)- gravity 

We consider a Bianchi spacctimcs whose homogeneous hypersurfaces have isotropic 
3-curvaturc ^Rab ~ ^{^R)hab- These spacetimes include the Bianchi models which, 
via the dissipation of the shear anisotropy aab, can reach a FLRW limit. Spatial 
homogeneity implies that the spatial gradients will vanish and that iia = = w. Thus 
the trace free Gauss-Codazzi equation (|16(l becomes 

^T(afc) + e (Tab = [^Tab + [h^ht - Ihabh'") S,d] , (18) 

and 5*06 can be split as follows 

Sab = f" (kuaUb - RVbUa) + f'R^UaUl,, (19) 

Sabu''u' = f"R + ,rR^, (20) 

Sabh'''' = -f'RQ, (21) 

S = - f"{R + m)- f"'R^. (22) 

Substituting these components into the Gauss-Codazzi equation p8|) gives 

O'(at) + Q CTab = f T^ab - f"Rcr{ab) ■ (23) 
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In the case of a perfect matter fluid. tTq^ = 0, so that the equation above becomes 

d f R 

O'iab) + 9 (Tafc = " ("^ ^at) = -JJ-<yab- (24) 

On integration this yields 

aab = f''^-^aba-\ *a6 = 0, (25) 
which in turn implies 

0-2 = /' -2^f2a-6^ ^2^0. (26) 

In the case of f{R) = i?, equation (|26|l gives the standard GR solution (see [11] and 
references there in) whose behaviour can be summarised as follows: 

(T'2 ^ oo as a — > 0, cr^ ^ as a —> oo. 

This behaviour is modified in /(i?)-theories of gravity (see [9,23]), because R (and 
therefore /') is a function of (see (|3Uf) below) and therefore H2t)|) is implicit. 
In particular, the dissipation of the shear in Bianchi I spacetimes is slower in 
quadratic gravity than in GR [9]. However, this result was obtained by solving the 
evolution equations under the assumption that the scale factor also has a power- 
law evolution. Although this is desirable it may not necessarily be true since no 
analytical cosmological solution could be obtained in [9] . A more general approach to 
this problem is to make use of the theory of dynamical systems (see [6] and references 
therein). In the following we will apply this technique to i?" gravity in order to 
investigate further the behaviour of the shear in this framework. 

3.2. R"- gravity 

We begin by specialising all the evolution equations above to the case of f{R) = -R". 
The Raychaudhuri equation H15|) is now 



e + ie^ + 2.^--i?-(.-i)-e + -^ = o, (27) 

and the trace free Gauss-Codazzi equation (|23|1 for LRS spacetimes is given by 



e + {n-l)^]a. (28) 



R 

The Friedmann equation can be found from (|17|l 

ie2-a2 + (n-l)|e-^^^i?--f^ = 0. (29) 
3 ^ 'R 2n nR'^-^ ^ ' 

In general, the substitution of the Friedmann equation (|29|l into the Raychaudhuri 
equation (|27() yields 

i? = 2e + |e2 + 2cr2. (30) 

Note that, in this relation the energy density does not appear explicitly, but is however 
still contained implicitly in the variables on the right hand side. 

In this paper we will assume standard matter behaves like a perfect fluid with 
barotropic pressure p = w[i. The conservation equation in this case is 

/( = -(l + w)^e. (31) 
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In order to convert the equations above into a system of autonomous first order 
differential equations, we define the following set of expansion normalised variables |; 



--l^i^-^)^ (32) 



3R 3fj. 



whose equations are 

E' = 2 (^-2 + 21] - -2x + z^ E, 

x' ^y(2 + x) — (2 + nx) - 2x - 2x'^ + xz + (1 - 3w)z + 2xT,, (33) 

71 — f 

y' = [(3 - 2n)x - 2y + 2{n - l)z + 4(n - 1)S + 2{n - 1)] , 

n — I 

2z - (1 + 3w) -3x — + 41] 

n — I 

where primes denote derivatives with respect to a new time variable r = In a and the 
dynamical variables are constrained by 

l-T. + x- y- z = 0. (34) 
4. Dynamics of the vacuum case 

We first consider the vacuum case (/i = 0). In this case the set of dynamical equations 
are given by 

E' = 2 (^-2 + 2E ^ - 2x^ E, 

x' ^y(2 + x) — (2 + nx) - 2x - 2x^ + 2xT,, (35) 

71—1 

y' = [(3 - 27i)a; - 2y + 4(n - 1)E + 2(n - 1)1 , 

11 — 1 

together with the constraint equation 

l-T. + x-y^O. (36) 

Fixed points and solutions 

The two most useful variables are E and y since they respectively represent a measure 
of the expansion normalised shear and the expansion normalised Ricci curvature and 
hence allow us to investigate how the shear is modified by the curvature. We can 
therefore simplify the system (|35|) . by making use of the constraint (|36|l . which allow 
us to write the ec^uation for a; as a combination of the two variables E and y: 

y' = [(211 - 1)E - (211 - 1)7; + (4n - 5)] , (37) 

n — 1 

J It is important to note that this choice of variables will exclude GR, i.c the case of n = 1. See [6] 
for the dynamical systems analysis of the corresponding cosmologies in GR. 
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Table 1. The fixed points and eigenvalues for _R"-gravity in a LRS Bianchi I 
vacuum model. 



Fixed points (S,y) Eigenvalues 



Points (0,m) 
Line £1 (i;*,0) 



2(5-4n) (5-4n) 
n— 1 ' n— 1 



(4n-5) , (2n-l), 



which together with the constraint H^{6|l represents our new system. Setting S' = 
and y' = we obtain one isotropic fixed point A : ^0, 2n-i ^ ^^'^ ^ -'^^^'^ '^^ fixed points 
£1 : (S*, 0) where S* > (S* < would imply imaginary shear) §. The point = 
on £1 represents another isotropic fixed point that merges with A when n = 5/4. 

The fixed points may be used to find exact solutions for the Bianchi I models. 
We substitute the definitions H32|l into H3U|) to obtain 

e=(^.,-E.-2)^, (38) 

where (S^jj/^) represents the coordinates of the fixed points. Given that n 7^ 1 and 
nyi — {n ~ I) (E^ + 2) 7^ 0, this equation can be integrated to give 

a^aoit-toT, where a ^ (2 + j:, - ^y,^ \ (39) 
In the case of the fixed point A we have 

(l-n)(2w-l) 

a = ao{t-to) ("-2) ^ (40) 

which is the same as the solution found in [7] . 
For the fixed line £1 we have 

a^ao{t-to)^ , (41) 

but direct substitution into the cosmological equations reveals that this solution is 
only valid for n > 1. For n < 1 the fixed points on £1 are non physical because the 
field equations do not hold there. We have summarised these solution in Table |21 

The analysis above would be incomplete without determining the fixed points at 
infinity. In order for us to compactify the phase space, we transform our coordinates 
(E, y) to polar coordinates 

E = fcos^, y = fsini^ (42) 

and set f = j^. Now since E > 0, we will only consider half of the phase space, i.e. 
— 7r/2 < 4> < tt/2. In the limit r ^ 1 (f — > 00), equations H37|l take on the form 

(2n — 1) 

— ' — s — cos 30 — 5 sin — sin 30 ] , (43) 



4(n- 1) ' 

I {2n — 1) [cos (f) — cos 30 + sin + sin 30 
= (n-l)(l-r) 
S* are the coordinates on the S-axis. 



(44) 
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Table 2. The solutions of the scale factor and shear evolution for i?"-gravity in 
a LRS Bianchi I vacuum model. 



Scale factor Shear 



(l-rt)(2K-l) 

Points a = ao(t-to) (""^^ cr = 

1 

Line Ci a = ao{t — to) 2+s, ^ (only valid for n > 1) a = croa"'-^'''^*-' 



Table 3. Asymptotic fixed points, 0-coordinatcs and solutions for ij" -gravity in 
a LRS Bianchi I vacuum model. 



Point (j> Scale factor Shear 



^oo T-Too = ^ln|i-<o| + C cr = (Toa-(2+Se) 



Boo 


TT 

2 


|t - Tool = 


Cl ± Co 


n-l 
2ri-l 




to) 


71 -1 

2n-l 


(7 = 


c 


3tt 
1 


|t - Tool = 


Cl±Co 


n-1 
2n-l 




to) 


n-1 


(7 = 




4 


|t - Tool = 


Cl±C2{ 








(7 = (70 



Since H43|) does not depend on r we can find the fixed points by making use of 
H44|l only. Setting i^' = we obtain four fixed points which are listed in Table 3 with 
their corresponding solutions. 

The form of the scale factor can be determined from the fixed points by integrating 
^ to find [8,24] 

r-l=g((/),)(T-Too), (45) 

where g{4>i) represents the right hand side of (|43|l and r Too as r — !■ 1. The evolution 
equation H38|l is then transformed into polar coordinates 

e' r / n 2(1 -r)\ 

— ' sin — cos (jji 1 , (46) 



e 1 - r \n - 1 
which in the limit r I take on the form 

%^-rr-\7^ ( ^^-sin(^j -cos(/), ) . (47) 

The equations above were all given in terms of the new time variable t by using the 
relation 0' = 30/9. Integrating (|47|l yields the solution 

|t - Tool = [Cl ± Co \hi^^)\ [t ~ to)] , (48) 

where 

'*(<^») = -TTT ( — — r - cos(/., ) + 1. (49) 

g{(j)i) \n-l ) 
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Solutions at infinity can now be obtained by directly substituting the fixed points into 
H48|l . so for points B^o and Coo we find 



n — 1 



{t - to) 



2n-l 

n-1 



(50) 



2n-l 
and for point 

\T-T^\ = [Cl±C2{t-to)f . (51) 

The solution at point Aca can not be determined with this method since the limit 
approaches our fixed line Ci so that l|48|l yields an indefinite solution. However, 
defining the two new variables: 5* = InS and Y = hiy, the system (|H7|I can be written 
as 



S' = - 2 



n — 1 
2n- 1 



n-1 



4n - 5 
n — 1 



(52) 



Point Aca corresponds to ?/ — > as S — > oo, so in its neighborhood the system ^ 
reduces to 

271- 1 



S" = and Y' 
which has the solution 



n-1 



S = Sc = constant and Y 



2n- 1 
n-l 



e^" (r-Too). 



(53) 



(54) 



The form of the scale factor for can then be found in the same way as the previous 
points. For y ^ as S — - ^-i'— ^ — ^ - ^ ^ 

find the solution as 



oo, H38|l takes the form 8 = — where 



-^ln|i-io| + C, 



e'^". We 
(55) 



where C is a constant of integration. 



4-2. Stability of the fixed points 

The stability of the fixed point may be determined by linearising the system of equation 
H37|l . This can be done by perturbing E and y around the fixed points (Si^yi) via 
E = Ei + (5E and y = yt + Sy. The corresponding eigenvalues of the linearised system 
are given in Tabled The fixed point A is an unstable node (repeller) for values of n 
in the range 1 < n < 5/4. For all other values of n it is a stable node (attractor). 

The fixed points on line £i all contain at least one zero eigenvalue and therefore 
we will have to study the effect of small perturbations around the line. We find that 
they have the following solutions 

ST,=^-e"\ Sy^Ce"^ (56) 

V 

where C is a constant of integration and 

(2n-l) , (4n-5) 2(2n-l) 
[n — 1) [n — 1) (n — 1) 
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Table 4. Stability of the fixed points for _R"-gravity in a LRS Bianchi 1 vacuum 
model. We use the term 'attractor' to denote a stable node and 'repeller' refers 
to an unstable node. 



n<l/2 l/2<n<l 1< 7i< 5/4 n > 5/4 



Point A attractor attractor repcUcr attractor 



Line Ci 

= 1 repeller repeller repeller repeller 

< < 2n-^ repeller repeller attractor repeller 

^* ^ 2n-i repeller attractor repeller repeller 



Table 5. Stability of the asymptotic fixed points for R^-gravity in a LRS Bianchi 
1 vacuum model. We use the term 'saddle' to denote a saddle node. 



Point n < 1/2 1/2 < n < 1 



n > 1 



r 



repeller 
saddle 
saddle 

attractor 



saddle 
saddle 
saddle 
repeller 



repeller 
saddle 
saddle 

attractor 



In order for the fixed points on line Ci to be stable node, wc must have < 0. The 
fixed point is an unstable node when 77 > 0. Over the interval < < , we will 
have stable nodes for 1 < n < 5/4 and unstable nodes for 1/2 < n < 1. The remainder 
of the points > 2n-i ' ^^^^ stable nodes for 1/2 < n < 1 and unstable nodes for 
1 < n < 5/4. When n < 1/2 and n > 5/4, the fixed points are always unstable nodes. 
Wc also note that for = 1, 77 = 6 and is therefore always an unstable node. The 
stability of all the fixed points is given in Table 4. 

A similar analysis may be performed for the fixed points at infinity. Wc only 
need to perturb the angular variable (f> around the fixed points 0^ via (j) = + S4>. 
The fixed points will be stable if r' > and the eigenvalue A < for the linearised 
equation S(p' = A 5(f), in the limit of f 00. When both conditions are satisfied the 
point is an stable node, if only one is satisfied it is a saddle and when neither holds 
it is an unstable node. Substituting the expression above into l|44|) and linearising as 
before, yields 

S(j)' « w ^^"i x/^^ — r [ - sin 0i + 3 sin 30i + cos (j>i + 3 cos 30i 1 Scj). (58) 
4(n — 1)(1 — r) 

The stability of the fixed points are summarised in Table 5. We see that only point 
I? 00 have stable nodes for n < 1/2 and n > 1. Points Boc and Coo are always saddle 
points and Aoc is a saddle when 1/2 < n < 1 but is otherwise an unstable node. 
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4-3. Evolution of the shear 

In the previous section we found two isotropic points; the fixed point A and one point 
on the fixed fine at = 0. The remaining fixed points aU have non-vanishing shear. 

The trace free Gauss Codazzi equation (|28|l can in general (i.e. for all points in 
the phase space) be represented in terms of the dynamical variables 1)32(1 as 

^ = ~i(2 + I] + y)e. (59) 

From the equation above it is clear that the shear evolution for all points in the phase 
space that lie on the line ?/ = 1 — S, is the same as in the case of GR. The shear will 
dissipate faster than in GR when a/a < —8, that is all points that lie in the region 
y < 1 — E. We will call this the fast shear dissipation (FSD) regime. When a/a > —Q 
and hence for all points in the region y > 1 — E, the shear will dissipate slower than 
in GR. This will be called the slow shear dissipation (SSD) regime. 

The fixed points on Ci for which E* > all have non-vanishing shear. For these 
points H5H|I has the form 

^ = -i(2 + E,)e = -(2 + E,)(^^^, (60) 

which may be integrated to give 

a = aoa-^'+^'^ = aoa-^'+'''\t - to)-\ (61) 

where we made use of (|41ll . We note that the final solution of the shear for these fixed 
points (|61|l . does not depend on the parameter n. This is to be expected since both 
the coordinates of the points on Ci and equation (|59() are independent of n. 

The only other fixed point with non-vanishing shear is I?oo which corresponds to 
y —> —oo as E ^ oo. In this limit ((59|l yields a/a = which implies that a = ao (i.e. 
constant shear). 

We first consider values of the parameter for which n < 1/2 (see Figure 1). If 
the initial conditions of the universe lie in the region y > (negative Ricci scalar), 
the orbits will always approach the isotropic fixed point A. The shear will dissipate 
slower than in GR for almost all the orbits in this region, apart from the ones below 
the dotted line y = 1 — S, which make the transition from the SSD region to the FSD 
region. When the initial conditions lie in the region y < (positive Ricci scalar), the 
orbits will approach the fixed point I?oo which has constant shear. 

The case 1/2 < n < 1 is illustrated in Figures 2 and 3. If the initial conditions 
are such that they lie in the shaded area {y < 1 — E and y < 0), then the evolution 
will always be in the FSD regime and will approach the isotropic solution of the point 
A. Instead, the unshaded area y < and y > 1 — E, is divided into two regions; the 
first one is located below the dash-dotted line and the second above the dash-dotted 
line. For initial conditions that lie in the first region, the orbits make a transition from 
the SSD region to the FSD region where they approach the point A. When the initial 
conditions lie in the second region, the orbits will always lie in the SSD region and 
approach Ci. When the initial conditions lie in ?/ < 1 — E and y > 0, the universe will 
evolve from the FSD regime to the SSD regime, in which the evolution will approach 
the stable solutions on Ci. In the remaining area where y > 1 — E and y > 0, the 
shear will always dissipate slower than in GR. 

We next consider 1 < n < 5/4 which is illustrated in Figure 4. If the initial 
conditions of the universe are such that they lie in the shaded area, then the shear 
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will always dissipate faster than in the case of GR. When the initial conditions lie 
in the region y > 1 — H and y > 0, the universe will evolve from the SSD regime to 
the FSD regime and the evolution will approach the stable solutions on Ci. If the 
initial conditions lie in the unshaded area of y < 0, the orbits will approach the fixed 
point Poo- For initial conditions that lie in the region y < 1 — E, there will be a 
transition from the FSD region to the SSD region. For all initial conditions that lie 
in the region y < 1 — S, the orbits always lie in the SSD region. We can see that 
for this range of n, the fixed point A acts as a past attractor. This is an interesting 
feature since an isotropic past attractor implies that unlike GR, where the generic 
cosmological singularity is anisotropic, we have initial conditions which corresponds 
to a FLRW spacetime. This feature was also found in the braneworld scenarios where 
it was shown that homogeneous and anisotropic braneworld models (and some simple 
inhomogeneous models) have FLRW past attractors (see e.g. [18-21]). This means 
that although inflation is still required to produce the fluctuations observed in the 
cosmic microwave background (CMB), there is no need for special initial conditions 
for inflation to begin [25]. In the range 1 < n < 5/4, we can obtain models whose 
evolution starts at the isotropic point A and then either evolves toward the fixed points 
(Eh., 0) on line £i or towards the point Doo- The orbits that approach £i will always 
lie in the FSD region. The orbits which approach Poo will make a transition from the 
FSD region to the SSD region. An interesting set of orbits are the ones that approach 
the fixed points on £i for which (cr/i?)* << 1. These cosmic histories represent an 
universe that is initially isotropic and then develops shear anisotropics which approach 
a constant value that can be chosen to be comparable with the expansion normalised 
shear observed today ((cr/iJ)* < 10~^ [26-28]) ||. Inflation is therefore not required 
to explain the low degree of anisotropy observed in the CMB. Furthermore, we still 
require all other observational constraints to be satisfied. 

Finally, we consider values in the range n > 5/4 (see Figure 5). For all initial 
conditions that lie in the shaded region, the shear will always dissipate faster than 
in the case of GR; for y > the orbits will approach A and for y < approach the 
point "Doo- If initial conditions lie in the region y > 1 — T, and y > 0, the orbits will 
initially be in the SSD region, and then approach the isotropic solution A, which is 
in the FSD region. For all orbits in the region y > 1 — E and y < 0, the shear will 
always dissipate slower than in the case of GR. 

5. Dynamics of the matter case 

We will now consider the dynamics of LRS Bianchi I models in the presence of matter. 
As noted in [7], in HOTG there is a difference between vacuum and non- vacuum 
physics in the sense that not all higher order couplings are consistent in the presence 
of standard matter. This can be seen in the evolution equations (|27() and H29|l where 
the matter terms are coupled with a generic power of the curvature. Since the sign 
of the Ricci scalar is not fixed, these terms will not be defined for every real value 
of n. Thus, the inclusion of matter induces a natural constraint through the field 
equations on i?"- gravity and it is therefore necessary to express the results in terms 
of the allowed set of values of n. Following [7], we will work as if n is unconstrained, 
supposing that the intervals we devise are meant to represent the subset of allowed 

II Strictly speaking these orbits do not satisfy the CoUins and Hawking [29] definition for 
isotropisation, which require a/H to asymptotically approach zero. 




Figure 3. Close-up of the phase space of the vacuum LRS Bianchi model with 
1/2 < n < 1 around the line y = I — T,. 




0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3 



Figure 4. Phase space of the vacuum LRS Bianchi model with 1 < n < 5/4. 
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Figure 5. Phase space of the vacuum LRS Bianchi model with n > 5/4. 



values within these intervals. 

Similar to the vacuum case, we can reduce the system to the three variables 
S, y and z by making use of the constraint (|34|l : 

' '2n-V 



1 



y 



1 



y + z 
[{2n- 1)1]- (2n 



(2 - 3w) - z + S 



1)?; + z + (4n - 5)] 
3n - r 



(62) 



We note that when y = then y' = and when 2 = 0, 2' = 0. The two planes y = 
and 2 = therefore corresponds to two invariant submanifolds. When 2 = 0, the 
system H62|) reduces to system (I37II and one would be tempted to consider the plane 
z = as the vacuum invariant submanifold of the phase space, which would not be 
entirely correct. To illustrate this point, we write the energy density in terms of our 
expansion normalised variables (|32|l 

^i oc zy"-ie"". (63) 

From this relation it can be seen that when z = and y ^ the energy density is zero. 
However when y = and z 7^ the behaviour of /i does depend on the value of n. In 
this case the energy density is zero when 71 > 1 but is divergent when n < 1. When 
both y and z are equal to zero and n < 1, one can only determine the behaviour of fi 
by direct substitution into the cosmological equations. 



5.1. Fixed points and solutions 

Setting E' = 0, ?/' = and z' = we obtain three isotropic fixed points A, B, C and 
a line of fixed points £1 : (E*, 0, 0), where > (sec Table 0). When = we 
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have another isotropic fixed point which merges with A when n = 5/4 and with B 
when w — 2/3. This point will merge with C when n = 5/4 and w ~ 2/3. 
We again substitute the definitions H32|) into (|30|) to obtain 

e^(^,.-E.-2j-. (64) 

Under the condition that n =/= 1 and the terms inside the brackets are not equal to 
zero, this equation may be integrated to give the following solution 

a^aoit-toT, where a = (2 + - ;^?/,) \ (65) 

The point A and line Ci will all have the same solutions as in the vacuum case, since 
either y = or z = for these points. 

The behaviour of the scale factor for point B is 

a = flo (t - to)^^^ , (66) 

and since 2/2 = and Z2 7^ 0, the energy density is zero (only valid for n > 1). When 
n < 1, point B and the points on Ci are non physical since the energy density is 
divergent. 

For point C, the scale factor behaves as 

2n 

a^ao{t~to)^(^ , (67) 
while the energy density is 

M = Moi"'", (68) 

where 

fj.0 = (-l)"3-"22"-in"(l + w)-^'\An - 3(1 + w))"-^ 
X [271^(4 + 3w) - ri(13 + 9w) + 3(1 + w)] . 

This point thus represents a power-law regime which in the case of 71 > 0, yields an 
expanding solution with the energy density decreasing in time. In the case of n < 
we obtain a contracting solution with ^ increasing in time. In order for C to be a 
physical point, we require /i > and therefore fj-o > (see [7] for detailed analysis). 
Note than when n > |(1 + lu), this solution corresponds to accelerated expansion. 

We next study the behaviour of the system H62|) at infinity. The compactification 
of the phase space can be achieved by transforming to spherical coordinates 

E — fsin0cos(/), y — fsm6 sin(f>, z = fcos0, (69) 

and setting f = jrj., where < f < 00, < 9 < tt and since we are again only 
considering half of the phase space, —tt/2 < cf) < tt/2 . In the limit r — > 1 (f ^ 00), 
equations H62|) take on the form 

/ = ^ [ 2cos0{3- 2n+ (2n- l)cos20} + 2 cos^ 6'{(2?i - 1) cos2(/i - 1} 

+2 sin 6* cos'^6{{2n- 1) cos 20- l}{cos<j> + sincf)} 
+ {2n — 1) sin 6' {cos0 — cos 30 — 5 sine/) — sin 30} ] , (70) 
, sin26'{l — (2n — 1) cos 20} [cos 61 + sin6'{cos0 + sin0} ] 



4(n - 1)(1 - r) 
(2n — 1) sin 20 [cos 6 + sin 6'{cos + sin 0} 
2(?i- 1)(1 -r) 



(71) 
(72) 
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Table 6. Fixed points and solutions of the scale factor and matter density for 
LRS Bianchi I with matter. 

Point Fixed points (S, y, z) Scale factor Matter density 

(l-n)(2n-l) 

A (O, 1^, Oj a = ao{t~to) (""2) m-0 

B (0, 0, 2-3w) a = ao{t~to)^ /i = (n > 1) 

(„-i)[4n^--3(i+,.)] ^ a = ao(t-to)3(i+-) M = Mor 



C 



j--2n 



2^ 



Line£i (S,,0,0) a = aa (t - to)^+^' M = (n > 1) 



Table 7. Eigenvalues and shear solutions for LRS Bianchi I with matter. 



Point Eigenvalues 



Shear 



A 
B 

C 



2(5-4rt) 5-4n n{13+9w)-2n^{4+3w)-3{l+iu) 
n-1 ' n-1 ' l-3n+2ra^ 
4w-3(l+ii)) 
n-1 



-2 + 3w, -4 + 6u;, 

3((27i-l)-m-l) Pi(n,tu)-^P2(",t") Piin,w) + y' P2(n,w) 
n ' 4n{n— 1) ' 4n(n— 1) 

Pi(n,u)) =3(l + ?«) + 3n((2n-3)w-l) 
P2(", w) = (n - 1) [4n3(8 + 3w)2 - 471^(152 + 3w(55 + 18w)) 
+3?7.(1 + w)(139 + 87w) - 81(1 + wf] 



cr = 
cr = 

cr = 



Line £i 



Now since ()70|l does not depend on r we can find the fixed points by just making 
use of (|71|) and H72|) . Setting 6*' ~ and 0' = we obtain the fixed points which are 
listed in Table with their corresponding solutions. 

The form of the scale factor can be determined from the fixed points as in the 
vacuum case. We integrate ((70)1 to find 

r-l^G{e,,cj),){T~T^), (73) 

where G{di,<j>i) represents the right hand side of H70() and r ^ Too as r ^ 1. We 
transform the evolution equation (|64|l to polar coordinates and write them in terms 
of the time parameter r: 

e' r sine, ( n 2(1 -r)\ 

sm (pi - cos <pi — - , (74) 



(75) 



e i-r V"-i 

which in the limit r ^ 1 takes the form 
e' -sin6',- 
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Table 8. Coordinates and solutions of the asymptotic fixed points for LRS 
Bianchi with matter. For the fixed line = arctan[— l/(cos cjii + sincpi)] and 
for the double point = (1/2) arccos[l/(2n — 1)]. 



Point 


(OA) 


Scale factor 










Shear 


A 

Boo 
(^oo 


(0,0) 

(77,0) 

(f.o) 


\t-Too\ = [Ci±Co(t-to)] 
\t-Too\ = [Ci±Co(t-to)] 
r-Too = ^ln|t-to| +C 


2TI-1 


cr = 
cr = 


Voo 


/ 7r 7r \ 
V 2 ' 2 / 


It-TooI = 


Ci±Co 


n-1 
2n-l 


- h) 


n-1 


(7 = 


f 


(IT 37T \ 

V 2 ' 2 > 


|t - Tool = 


Ci±Co 


n-1 
2n-l 


(t-ta) 


271-1 

n-1 


cr = 


oo 




It-TooI = 


Cl±C2( 








cr = (Jo 


Line 
















r 

*^oo 


{o,4>) 




Ci±C2(t-to)]' 




(7 = CTo 



Integrating the equation above then yields 



[Ci±Co|i/(0.,4)l(i-to)]^(''-^') 



where 



snif 



■ smf 



C0S( 



1. 



(76) 



(77) 



The solutions is then obtained by directly substituting the coordinates of the fixed 
points into H76() . For points Aoo and Boo we have the solutions 

|r-roo| = [Ci±Co(i-M], (78) 
and for point I?oo and £oo we have 



n- 1 



2n- 1 



{t - to) 



(79) 



In addition to the ordinary fixed points, we have a line of fixed points for which 9 = 
arctan[— l/(cos(/ij; + sini/)^)], and a double fixed point a,t (p = (1/2) arccos[l/(2n — 1)]. 
Equation (|76|l gives the solutions for the fixed points on the line: 

\T-Too\^[Ci±C2it-to)f . (80) 

As in the vacuum case, Coo can not be determined with this method since the 
limit approaches the fixed line Ci. In order to analyse this point, we can define three 
new variables: S* = InS, Y = lny and Z — \nz. The system ()62f) can then be written 
as 

'2ri- 1' 
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Table 9. Stability of the fixed point A. for LRS Bianchi I with matter. The 
parameters are N± = ^(13 ± V73), P± = ^{A ± VE) and Q± = ^(11 ± V37). 







n < N- 


N- < n < P- 


P < n < O 


O < n < 1/2 


w 


= 


cittractor 


saddlc 


saddle 


saddle 


w 


= 1/3 


attractor 


attractor 


saddle 


saddle 


w 


= 1 


attractor 


attractor 


attractor 


saddle 






1/2 < n < 1 


Kn <Q+ 


Q+ <n< 5/4 


5/4 < n< P+ 


w 


= 


attractor 


repeller 


repeller 


saddle 


w 


= 1/3 


attractor 


repeller 


repeller 


saddle 


w 


= 1 


attractor 


repeller 


saddle 


attractor 






P+ <n<N+ 


n> N+ 






w 


= 


saddle 


attractor 






w 


= 1/3 


attractor 


attractor 






w 


= 1 


attractor 


attractor 







The point Coo corresponds to y — > and z — > as S — > oo so that the system ()81|l 
reduce to 

S' = 0, Y'^ (^TTt) ^' = (82) 

which has the foUowing solution 

S = Sc, y= (^^^^^e^^(r-Too), and Z = e^^(r-Too),(83) 

where Sc is a constant. The form of the scale factor for Coo can then be found in the 
same way as the previous points. For y ^ and z — ^ as S ^ oo, (|64|l takes the 
form = — where Sc = e'^''. We find the solution as 

T-Too = ^\ll\t-to\+C, (84) 

where C is a constant of integration. 
5. 2. Stability of the fixed points 

We next check the stability of the fixed point by linearising the system of equation 
Ht)2|l . The eigenvalues of the linearised system are given in Table 7. 

For the stability analysis, we consider three cases: dust w = 0, radiation w — 1/3 
and stiff matter w = 1. The stability of the fixed points A and B are summarised 
in Tables ini and ^1 respectively. Their behaviour is similar to the flat (fc = 0) points 
(C and J^) considered in [7]. The stability analysis for the fixed point C cannot be 
performed in an exact way. The eigenvalues are complex in the following ranges: When 
w; = 0, for 0.31 <n< 0.71 and 1 < n < 1.32; when w = 1/3, for 0.35 < n < 1.28; and 
when w — 1, for 0.37 < n < 1 and 1.224 < n < 1.47. The results are given in Table 

m 
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n<3/4 3/4<7i<l 1< n < 3/2 n > 3/2 



u) = saddle attractor saddle saddle 

w = 1/3 saddle saddle saddle saddle 

w — 1 rcpellcr repeller saddle repcUer 



Table 11. Stabihty of the fixed point C for LRS Bianchi I with matter. We use 
the term 'spiral' to denote a pure attractive spiral and an 'anti-spiral' refers to a 
pure repulsive spiral. 







< n < 0.31 


0-31 < n < 0.35 


0-35 < n < 0.37 


0-37 < n < 0.71 


0.71 <n< 0.75 


w — 





attractor 


spiral 


spiral 


spiral 


attractor 


w — 


1/3 


attractor 


attractor 


spiral 


spiral 


spiral 


w = 


1 


saddle 


saddle 


saddle 


spiral 


spiral 






0.75 < n < 1 


1 < n < 1.220 


1.220 < n < 1.224 


1.224 <n< 1.28 


1.28 < n < 1.29 


w — 





saddle 


spiral 


spiral 


spiral 


spiral 


w — 


1/3 


spiral 


spiral 


spiral 


spiral 


attractor 


w — 


1 


spiral 


saddle 


repeller 


anti-spiral 


anti-spiral 






1.29 < n < 1.32 


1.32 < n < 1.34 


1.34 < n < 1.47 


1.47 < n < 1.5 


n > 1.5 


w — 





spiral 


attractor 


saddle 


saddle 


saddle 


w — 


1/3 


saddle 


saddle 


saddle 


saddle 


saddle 


w — 


1 


anti-spiral 


anti-spiral 


anti-spiral 


repeller 


saddle 



As in the vacuum case we find that the fixed points on the line Ci have zero 
eigenvalues. We therefore study the perturbations around the fixed line which lead to 
the following solutions 

where Co and Ci arc constants of integration and 

(2n - 1) Un - 5) 
(n-1) [n-l) 

In order for the fixed points on line Ci to be stable nodes, we must have r] < and 
K < 0. When 77 > and k < or 77 < and k > we have a saddle and when 77 > 
and K > it is an unstable node. The results have been summarised in Table [T^ 

A similar analysis can be performed for the fixed points at infinity. Wc can check 
the stability of the fixed point by linearising the system of equation H70|l - (|72|l . The 
eigenvalues of the linearised system are given in Tabled) The stability of the fixed 
points Aqo , Boo J Coo J T^oc and £00 can then be found straightforwardly as in the vacuum 
case (see Table Points Aoo and Boo are always saddle points. The point Coo lie 
on the fixed line £1 and is an unstable node for 71 < 1/2 and 71 > 1 and a saddle 
for 1/2 < 71 < 1. Point 2? 00 is an unstable node for 77 < and 71 > 1, a saddle for 
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Table 12. Stability of the fixed line Ci for LRS Bianchi I with matter. The 
coordinates for all these points arc (E*,0, 0) and therefore we only state the 
range of values on the fixed line in terms of their S-coordinates. 



w = 0, 1/3 Attractors Saddles Repellers 



n < 1/2 
1/2 < n< 1 
l<n<2 
71 > 2 



none 
none 
none 
none 



2n 

< E* < 



> 

< E. < 1^ 

2n— 1 

E. > 



w — 1 



n < 1/2 
1/2 < n< 1 



none 
none 



1< 7i< 5/4 < E, < 1^ 
n > 5/4 none 



< E, < 1 

< E. < 1 
^ < E. < 1 



< E, < 1 



E, > 1 
K E* < 



E. > 1 
E. > 1 



Table 13. Eigenvalues and value of r' of the ordinary asymptotic fixed points 
for LRS Bianchi I with matter. 



Point {9, (p) Eigenvalues r' 





(0,0) 


1 2ri-l 
' ri-1 




-1 


Boo 


KO) 


1 2n-l 
n-1 




1 


c 


(f'O) 


1 2n-l 
^' n-1 











/ 7r TT \ 

11 


n 271 

n— 1 ' n- 


-1 
-1 


2n-l 
n-1 




iTX 37r \ 
V 2 ' 2 ^ 


271-1 

n-1 ' 


n 


n 

-1 


2n-l 
n-1 



< n < 1/2 and a stable node for 1/2 < n < 1. Point is a stable node for n < 
and n > 1, a saddle for < < 1/2 and an unstable node for 1/2 < n < 1. 
The eigenvalues of the fixed line £00 are given by 



[Ai,A2] = 



Pi(n,0)- Pi(n,(/)) + VP2(n,0) 



4(7i- l)(2 + sin20)3/2' 4(n - 1)(2 + sin 2(/))3/2 



(87) 



where 



Pi (71,, 

P2(«,< 



(3?i — 4) cos (p + n cos 30 — {3n + 1) sin + (rt — 1) sin3(/), 

2(2 + sin 20)3/2 ^^27i2 - 2?^ + 1)(1 + sin 20) 

+ (2ri- l)(-cos20- 2sin40+ (2?7, - 1) sin 60)] . 

The stability can then be found in a similar fashion as the ordinary asymptotic fixed 
points. The eigenvalues in this case is dependent on two variables, n and 0, which 
makes it difficult to express the results in a table. We have therefore summarised these 
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Table 14. Stability of the ordinary asymptotic fixed points for LRS Bianchi I 
with matter. Results are independent of w. 



Point 


n < 


< n < 1/2 


1/2 < n < 1 


n> 1 




saddle 


saddle 


saddle 


saddle 


Boo 


saddle 


saddle 


saddle 


saddle 


^oo 


repeller 


repeller 


saddle 


repeller 




attractor 


saddle 


saddle 


attractor 




repeller 


saddle 


attractor 


repeller 




-1.5 -1 -0.5 0.5 1 1.5 

Figure 6. Stability of the line of fixed points Coo for LRS Bianchi I with matter. 
Results are independent of w. The black dot represents the fixed point at <^ = 0.4 
for a model with n = 1.4. 



results in a diagram (see Figure 6) %. The stability of any fixed point on the line Coo 
for a given value of n can be read from this diagram. For example the black dot in 
Figure 6 represents the fixed point at = 0.4 for a model with n = 1.4. It lies within 
a region that classify it as an attractor. 

1 This diagram was found by plotting r = Ai + A2, A = A1A2 and — 4 A and using the definitions 
for stability to classify the regions [30] . 
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5. 3. Evolution of the shear 

The trace free Gauss Codazzi equation (|28|) can in general (i.e. for all points in the 
phase space) be represented in terms of the dynamical variables (|32ll as 

- = -i(2 + S + y + z)e. (88) 

The shear evolves at the same rate as in GR when a /a = —6. which holds for all 
values on the plane 1 = S + y + z. The shear will dissipate faster than in GR when 
a /a < —6, that is all points that he in the region 1 > T, + y + z. We may again call 
this the fast shear dissipation (FSD) regime. The shear will dissipate slower than in 
GR when a /a > —8 and hence all points in the region 1 < E + y + z. This will be 
called the slow shear dissipation (SSD) regime. Analysing the three dimensional phase 
space for this system is more difficult than the two dimensional spaces considered in 
the vacuum case since it is harder to visualise. 

All the finite fixed points {A, B and C) together with the point = on Ci, lie 
on the plane E = and are therefore isotropic. 

The evolution of the shear for the anisotropic fixed points on Ci can be obtained 
as in the vacuum case. For these points H88|) take the form 

- = -i(2 + E,)e = -(2 + E,) f-V (89) 
a \a J 

which may be integrated to give 

where we made use of H39(l . This is the same solutions that where obtained in the 
vacuum case. This is to be expected since all these fixed points lie on the line Ci for 
which z ~ 0. 

The point Coo and the points on the line C^o (including J-oo) are the only 
asymptotic fixed points with non-vanishing shear. Point Coo lies on the line Ci and 
therefore has the solution 

The behaviour of the shear for the fixed points on Coo can be found like (|75|l . i.e by 
transforming (|88() into polar coordinates and taking the limit r ^ 1. The resulting 
equation then yields the solution a = ao — constant. 



6. Discussion and Conclusions 



We have derived the evolution equations of the shear for Bianchi I cosmologies with 
/(i?)- gravity. This general expression, allows us to consider the shear evolution for 
any function of the scalar curvature. However, because the shear depend non-lincarly 
on the Ricci scalar one can not determine how the dissipation of the shear anisotropy 
compares with the case in GR even if one chooses a specific form for f{R) (such as 
R + or i?"). One way of dealing with this problem is to make certain assumptions 
(for example the form of the evolution of the scale factor a [9]) to obtain a solution. 
A more general approach is to make use of the dynamical systems approach to study 
HOTG in these cosmologies since it provides both exact solutions and the global 
behaviour of the system. 

Our main aim in this paper was to see how the shear behaves in LRS Bianchi I 
cosmologies with i?"- gravity and whether these models isotropises at early and late 
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times. To achieve this goal we used the theory of dynamical systems to analyse the 
system of equations governing the evolution of this model with and without matter. 

The phase space for these models have a number of interesting features, in 
particular it contains one isotropic fixed point and a line of fixed points with non- 
vanishing shear. The isotropic fixed point is an attractor (stable node) for values 
of the parameter n in the ranges n < 1/2, 1/2 < n < 1 and n > 5/4. In the 
range 1 < n < 5/4 this point is a repeller (unstable node) and therefore may be 
seen as a past attractor. An isotropic past attractor implies that inflation can start 
without requiring special initials conditions. However, since we have attractors for 
(a/H)^ << 1 on £i, we may not need inflation since the shear anisotropy approaches 
a constant value which may be chosen as the expansion normalised shear observed 
today ((cr/i/)* < 10~® [26-28]), provided that other observational constraints such as 
nucleosynthesis are satisfied. 

We also found that the line y = 1 ~ T, separates the phase space into two part. 
For all points on this line, the shear dissipates at the same rate as in GR. In the region 
above the line the shear dissipates faster (FSD) than GR and in the region below the 
line, slower (SSD) than in GR. From Figures 1 — 5 we can see that there are a number 
of orbits which cross the dotted line. These are systems which initially lie in the FSD 
region and then make a transition to the SSD region and vice versa. An interesting 
feature of the vacuum case is that when the evolution of the universe reaches the stable 
solutions on £i, the shear will evolve according to a oc t~^ irrespective of the value 
of n. For values that lie in the range 1/2 < n < 1 and y > 0, one may have orbits 
that initially lie in the FSD region (see Figures 1 and 2) and then make a transition 
to the SSD region at late times. The opposite will happen for values that lie in the 
range 1 < n < 5/4 and y > 0- Initially the orbits lie in the SSD region (see Figure 3) 
and then make a transition to the FSD region. 

We observe the same kind of behaviour in the matter case where the phase space 
is however 3-dimensional, but is similarly divided into two regions, by the plane 
1 = T, + y + z. The space above the plane is the SSD region and below the FSD 
region. Similar argument to the vacuum case can be used here to investigate the 
orbits. When matter is included we do however only have stable fixed points on Ci 
for values of n in the range 1 < ti < 5/4. 

In conclusion we have shown that i?"- gravity modifies the dynamics of the shear 
in LRS Bianchi I cosmologies by altering the rate at which the shear dissipates. There 
are cases in which the shear always dissipate slower or faster than in GR, and there 
are ones which make transitions from first evolving faster and later slower (and vice 
versa) than in GR. 
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